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SUMMARY 


This paper describes the development and usage of an 
atmospheric Acoustic Source Bearing Estimation (ASBE) program. 
The basis for this program is a maximum-likelihood method of 
spectral analysis which was formulated by J. Capon for seismic 
array processing. 

Presented herein is the mathematical development of the 
Acoustic Source Analysis Technique (AS AT ) beari ng estimati on 
algorithm that is used in ASBE. Included in this report is the 
input and output of a test case that was used to validate the 
al gori thm. 
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SYMBOLS 


X ( i ) , Y ( i ) 


block averaged spectral estimate 
complex conjugate of B 

complex valued column matrix representing steering 
vectors 

conjugate transpose of E 

block averaged auto and cross spectrum 

complex conjugate of F 

equivalent to EE' 

hertz 

number of data points in each subblock of data 
number of data points being analyzed 
number of sensors 

number of blocks that L data points are separated 
into ( L=KN ) 

number of sensor pairings whose separation of 
distance is less than or equal to one-half the 
wavelength of the frequency being analyzed 

f requency-wavenumber power spectrum estimate 

ASAT power spectrum estimate 

complex valued square matrix representing the 
diagonally normal 1 zed-bl ock averaged estimate of the 
cross power spectral densities 

inverse of R 

complex values square matrix representing the cross 
power spectral densities obtained from the Fast 
Fourier Transform 

rectangular coordinates of sensor i 

speed of sound 

frequency being analyzed 

complex variable (7^1) 

time step between data points 

azimuthal look direction 
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INTRODUCTION 


Acoustic array processing entails the analysis of signals 
which are carried by propogating atmospheric sound waves. These 
waves impinge upon an array of sensors from which the received 
signals are recorded as pressure variations. Through array 
processing, the pressure values can yield signal characteristics 
such as power-frequency signature and source bearing. 

A major component of signal processing is the estimation of 
power spectral densities. Techniques which have been devised to 
estimate spectral values are divided into two major categories - 
window-based and model-based. The first category, window-based 
methods, are non-adaptive in the sense that all sets of data are 
treated alike. The most widely known and used window-based 
method is the Fast Fourier Transform (FFT). For a single sensor, 
the FFT will estimate the power-f requency signature of the signal 
being analyzed. From this signature, estimates can be made for 
the power spectral density and auto correlation values. When 
signatures are available from two sensors, estimates can be made 
for the cross power spectral density and cross correlation 
values. When a signal signature is required as output from the 
entire sensor array, the second category of spectral estimation, 
model-based methods, must be utilized. 

Unlike window-based methods, model-based methods of spectral 
estimation are considered adaptive since they are data 
dependent. For each set of data being analyzed, the 



characteristics of the model-based estimator are automatically 
adjusted. This adaptive capability allows much higher spectral 
resolution than the FFT technique. Since source bearing is of 
critical importance in acoustic analysis, the model-based method 
that was investigated is the maximum likelihood method (MLM). As 
documented by Capon*, this Maximum Likelihood Method provides a 
high-resolution estimate for frequency-wavenumber spectrum. The 
MLM acts as a filter at a selected frequency which passes that 
frequency undistorted and rejects all others. For a particular 
frequency, the MLM spectral estimates are a function of azimuth 
and are power levels that are received by the sensor array as a 
whole. Therefore, the azimuth which results in the largest power 
level is an estimate of source bearing. 

Both the MLM and FFT methods of spectral estimation have 
unique qualities which are beneficial in acoustic array 
processing. This report will present the derivation and 
implementation of ASAT, a technique to estimate acoustic source 
bearing which utilizes the capabilities of each method. 
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TECHNICAL DISCUSSION 


When the FFT is used to estimate each sensors' signal 
signature the highest frequency in the spectrum is equal to one- 
half the data sampling rate. This cutoff frequency is referred 
to as the Nyquist frequency. Additionally, the FFT is symmetric 
in the sense that there are an equal number of positive and 
negative frequency components in the spectrum. The frequency 
resolution of the FFT is therefore determined by this symmetry 
and the Nyquist frequency. For example, if L datum is sampled at 
time intervals of At and L datum is analyzed by the FFT then: 

Nyquist frequency = 1/(2 At ) = L/2 

frequency resolution = ( 1/2 At ) / ( L/2 ) = (L/2)/(L/2) = 1Hz. 

In order to increase the frequency resolution, augmenting zeros 
can be added to the sampled data. Using the example given 
previously, the Nyquist frequency remains l/2At and if L zeros 
are added to the data then: 

frequency resolution - ( 1/2 At ) / ( 2L/2 ) = (L/2)/L * 1/2 Hz. 

Thus, the frequency resolution has been doubled. An additional 

benefit of this zero padding is that foldover aliasing has been 
reduced. This type of aliasing occurs when frequencies higher 
than the Nyquist frequency occur in the data. These higher 

frequency components cause power foldover and distort the 
spectrum at lower frequencies. 

The multi-channel processing technique developed by Capon 
also requires spectral estimates for each sensor. Unlike the 
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FFT, the MLM utilizes a block averaging method to estimate a 
signal's signature. In this method a sample of a sensor's 

recorded signal is subdivided into blocks which contain an equal 
number of datum. The number of blocks must be greater than or 
equal to the number of sensors. For an array containing M 
sensors, each having a signal length of L data points which is 
subdivided into N blocks of K points such that N 2. M and L = KN, 
each sensor's block averaged spectral estimate in the nth block 
at a preselected frequency f is: 



IK)’ 1 ' 2 


K 

Z K 
m = l 


i , m+ ( n - 1 


im2 nf At 

)K e 


1-1,..., M (1) 

n-1 , . . . , N . 


The frequency resolution of the block averaging method is lower 
than that of the FFT. Whereas the FFT could have a resolution of 
1 Hz or less, this method would result in a resolution of only 
1 / ( K At ) . For example, if L = l/At = 2000, M=9 and N = 10, then K = 200 
and the block averaged frequency resolution would be 2000/200 = 
10Hz. 

Once spectral values have been estimated for each block in 
each sensor, the MLM requires that auto and cross spectrum be 
calculated by: 



N 

E 

n*l 




i , 1 = 1 , . . . M 


( 2 ) 


where B* indicates complex conjugate. A diagonal normalization 
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is then performed by dividing by These spectrum 

values are then arranged in a square matrix as: 



F 11 

F 12 

F 13 

• 

F 1M 

* 

F 21 

CM 

CM 

U. 

F 23 
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F 31 

F 32 

F 33 

• 

• 


F 41 

• 

• 

• 
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f mi 

• 

• 

• 

f mm 


where the li indexed values below the diagonal are the complex 
conjugates of the il indexed values above the diagonal. Once 
this complex valued matrix has been arranged, its inverse must be 
calculated. Situations can exist where the matrix is singular 
and thus has no inverse. If there is not enough reasonable data 
in each block or if the signal is transient in nature, then the 
matrix will be singular. When these conditions arise, the 
frequency under study and possibly the entire database will not 
be able to be analyzed using the MLM technique. 

If the inverse of the spectrum matrix exists, then the high- 
resolution estimate for the frequency-wavenumber array spectrum 
i s : 

P = [E ' R" 1 E]“ 1 (3) 
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where P is the real, positive valued power estimation of the 
spectrum; E is a column matrix whose complex elements represent 
steering vectors for each sensor; E‘ is the conjugate transpose 
of E; and R _1 is the inverse of the complex valued square matrix 
representing the diagonally normalized block averaged estimation 
of the cross power spectral densities at frequency f. Each 
element of the matrix E is of the form: 

E(i) « exp[-j2nf (X(i )C0S8 + Y(i)SINB)/c] (4) 

where f is the frequency, X ( i ) and Y(i) are the rectangular 
coordinates of sensor i, 8 is the azimuthal look direction, and c 
is the speed of sound for the local atmospheric conditions. 

Implementation of equation 3 is governed by the following 
assumpti ons : 

1. analysis is done on a frequency by frequency basis 

2. the output from a sensor is a wide-sense stationary 
discrete-time parameter random process with zero mean 

3. the output from the sensor array comprises a homogeneous 
random field 

Under these assumptions, equation 3 can be used to estimate the 
bearing of the source by determining which azimuth maximizes the 
power P at a selected frequency f. 

The preceeding discussion has presented a f requency-power 
spectral estimation technique, the rFT , and a multi-channel 
processing technique, the MLM . Each technique has qualities 
which are beneficial in acoustic array processing. Namely, the 


6 



FFT has high frequency resolution and anti-aliasing whereas the 
MLM will estimate source bearing. Spectral estimates obtained 
via the MLM were shown to be of low frequency resolution and 
could also give rise to a singular matrix. Therefore, it would 
be advantageous to incorporate the positive qualities of each 
technique into a unified model. The remainder of this paper will 
present the mathematical development of the ASAT algorithm and 
the implementation of this hybrid model of bearing estimation. 
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DERIVATION OF THE ASAT LOCATION METHOD 


The high resolution power estimate equation is: 


P * CE * R“ 1 E]“ 1 


since P is a scalar, the inverse of this equation is: 

1/P = E'R" 1 E 


multiply both sides by E'R: 

(E'R)/P = E'R“ 1 E{E'R) 


multiply both sides by E: 


(EE'R)/P = EE 1 R -1 (EE ' R) 


since P is a scalar: 


(EE’R) = PEE ' R" 1 (EE ' R) 


therefore: 


PEE'R -1 = identity matrix I 


( 5 ) 


(6) 


(7) 


( 8 ) 


(9) 


(10) 
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multiply both sides of equation 10 by R: 

PEE ' R _1 R = IR (11) 

therefore : 

P ( EE ' ) = R. (12) 

Equation 12 eliminates the need of computing the inverse of 
matrix R. Since the inverse is no longer required, matrix R need 
not be computed from the block a veraged-di agonal ly normalized 
method. Instead, the matrix can be obtained from the Fast 
Fourier Transform values. This then allows the use of the higher 
frequency resolution and anti-aliasing features of the FFT. 
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APPLICATION OF THE ASAT LOCATION METHOD 


Let the matrix G represent (EE*), matrix S represent the 
cross power spectral densities obtained from the Fast Fourier 
Transform, and P w represent the power estimate. Equation 12 can 
be rewritten as: 


P W G = S. (13) 

Since P w is a scalar value, every element of matrix S must equal 
every element of matrix G multiplied by P w . Therefore, P w must 
equal every element of S divided by each corresponding element of 
G: 


P w G(i,j) = S ( i , j ) => P w = $(i,j)/G(i,j). (14) 

In practice, the value of P w will vary since the S values 
are estimates obtained from a small segment of signal time 
history. To account for this variance of P w , an ensemble 
averaged value will be used: 

M M 

P w = i=l j=i (S(iJ)/G(i,j)) (15) 

1 * 5 ? 

The double summation over M sensors in equation 15 is 
constrained to i-j sensor combinations whose separation distance 
is less than or equal to one-half the wavelength of the frequency 
being analyzed. By applying this constraint, spatial aliasing is 


avoided. NSP is the number of sensor pairings which satisfy this 
distance constraint. 

Bearing estimation can be accomplished by using equation 15, 
the ASAT algorithm. 
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ASBE PROGRAM FLOW ANALYSIS 


The program which was developed to locate acoustic sources 
(Acoustic Source Bearing Estimation - ASBE) can be analyzed in 
six parts. 

These operational components are: 

1) read test site parameters from file MIKE 

2) read sensor output data from file TIFT 

3) convert sensor output data from dynes/sq-cm to Pascals 
and apply Hamming window 

4) compute cross power spectral densities by Fast Fourier 
Transform and compute ensemble averaged values for each 
frequency 

5) locate N peak ensemble averaged values where N is the 
degrees of freedom (number of sensors less 1) 

6) for each peak ensemble averaged value: 

a) compute power values at regular intervals of azimuth 

b) locate the peak values in the azimuthal power 
array, these values give the azimuths from which the 
source signals are assumed to radiate. 
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ASBE PROGRAM PARAMETER STATEMENTS 


Certain arrays 
parameter statements 

Vari able 
MAXNS 
NMIKES 
SPS 

DELAZ 


in the ASBE program are dimensioned by 
The user set parameter declarations are: 


Descri pti on 

integer number of sensor output values to read 
integer number of sensors 

real value for the sampling rate of the data 
(units: samples per second) 

real value defining the resolution to be used 
in the azimuth power calculation (units: degrees) 
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ASBE PROGRAM INPUT/OUTPUT 


Input to the ASBE program consists of the free formatted 
file MIKE and the unformatted file TIFT. 


File MIKE describes test conditions and contains: 


Record 

Variables 

Descri pti on 

1 

LABEL 

test name, maximum of 80 characters 

2 

NMK.AIRT 

integer value defining number 
of sensors used, real value 
defining test site air temperature 
in degrees Farenheit 

3 

PEAKI ,PEAKF 

real values defining the minimum 
and maximum values of the frequency 
range (hertz) to use in the peak 
search 

4 

X * Y 

coordinates for sensor #1 in 
feet 

5 

X,Y 
• » 
• • 

sensor #2 coordinates 

NMK+3 

• • 
X , Y 

sensor #NMK coordinates 


Note that the X,Y values defining each sensor position are 
in reference to a sensor which was selected as the origin and has 
a X-Y location of (0.,0.). 
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File TIFT contains the recorded pressure data and its 


structure is: 


Record 

Vari abl es 

Descri pti on 

1 

ISN.NC, (NAM( I ), 1 = 1, NC), 
(IU(I),I=1,NC),(IHD(J), 
J-1,8) 

integer test serial number 
integer number of channels 
recorded, 10-character name 
of each channel , 10- 
character units name for 
each channel, 80-character 
header name for the test 

2 

( CH ( I ) , I = 1 , NC ) 

• 

• 

C H { 1 ) : real-time value, 
CH ( 2 ) thru CH(NC): real 
pressure values for NMK 
sensors at time C H ( 1 ) 

MAXNS+1 

♦ 

(CH ( I ) , I = 1 , NC ) 

pressure values at ending 
time value 


The formatted output file from ASBE is named BEARING. It 
contains a list of program parameter values, an echo of input 
file MIKE, the calculated peak frequencies and associated power 
values, and the peak azimuth/power values for each peak 
frequency. An example of files MIKE and BEARING appear in 
figures 2, 3a, and 3b. 
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SAMPLE TEST CASE 


The acoustic source which was used to benchmark the ASBE 
program consisted of a stationary, gasoline fueled, auxiliary 
power unit (APU) equipped with a muffler. With respect to the 
sensor array origin (sensor #5), the APU was located at 120 

degrees azimuth and at a distance of 50 feet. The sensor array 
was comprised of nine B&K half-inch sensors arranged in three 
triangles with three sensors in each triangle. Figure 1 shows 

the test site geometry and sensor array setup. 

Sensor output data for this test was digitized at 2000 

samples per second for 30 seconds. For the benchmark analysis, 
this database was separated into 60 TIFT files. Each file 

contained 1000 samples covering a one-half second time block. 
Parameter statements setup for these 60 runs were: 


Parameter 

Val ue 

SPS 

2000.0 

MAXNS 

1000 

NMIKES 

9 

DELAZ 

1.0 


Each run was executed on a Control Data Corporation CYBER 
860 computer and required approximately 2.75 seconds of CPU time 
for each frequency analyzed. The input file MIKE used in these 
runs is shown in figure 2. 
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As a representati ve example of program execution, figures 3a 
and 3b show output file BEARING for the run covering the first 
one-half second of data. A histogram showing the typical 
variance of P w (obtained from equation 15) is presented in figure 
4. The power values used in this plot were calculated using a 
frequency of 140 hertz and a bearing of 121 degrees. 
Approximately 46 percent of the sensor pairings had power values 
that were within ± 0.6 percent of the mean. Figure 5 is a plot 
of the ensemble averaged spectrum magnitude. It is from this 
data that the peak frequencies used in this analysis are 
selected. The eight peaks selected for this run are highlighted 
by an asterisk. Polar power-bearing plots for each of the eight 
selected peak frequencies appear in figure 6. At each frequency 
the nominal direction for the maximum power is 120 degrees. It 
can be seen that sidelobes become more apparent at higher 
frequencies. This effect is caused by the sensor locations used 
in the test. As the frequency increases, more sensor pairings 
are excluded in the analysis since they have separation distances 
which are greater than one-half the wavelength of the frequency 
being analyzed. Thus, there is less data to use in the bearing 
estimation and this results in less discrimination of the 
sidelobes. In order for this effect to be minimized, the highest 
frequency would have to be known prior to the test so that the 
sensor array could be sized accordingly. 
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RECORD 

CONTENTS 

DESCRIPTION 


APU WITH MUFFLER 

label 


9 

56.0 

number of sensors, 

air temperature (degree Farenheit) 

50.0 

250.0 

frequency range 

(hertz ) 

-4.0 

0.0 

X-Y coordinates 
sensor 1 (feet) 

for 

4.0 

0.0 

X-Y coordinates 
sensor 2 (feet) 

for 

0.0 

-6.9 

X-Y coordinates 
sensor 3 (feet) 

for 

-2.0 

-3.4 

X-Y coordinates 
sensor 4 (feet) 

for 

0.0 

0.0 

X-Y coordinates 
sensor 5 (feet) 

for 

2.0 

-3.4 

X-Y coordinates 
sensor 6 (feet) 

for 

0.0 

-3.4 

X-Y coordinates 
sensor 7 (feet) 

for 

-1.0 

-1,7 

X-Y coordinates 
sensor 8 (feet) 

for 

1.0 

-1.7 

X-Y coordinates 
sensor 9 (feet) 

for 


Figure 2: 


Input File MIKE 
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***** RUN DIRECTIVES ***** 


PROGRAM PARAMETER STATEMENT SETS THE SAMPLES PER SECOND * 2000.0 

PROGRAM PARAMETER STATEMENT SETS THE AZIMUTH RESOLUTION (DEG) * 1.0 

PROGRAM PARAMETER STATEMENT SETS THE NUMBER OF DATA POINTS = 1000 

PROGRAM PARAMETER STATEMENT SETS THE NUMBER OF MICROPHONES * 9 

CASE LABEL : APU WITH MUFFLER 

NUMBER OF MICROPHONES : 9 

AIR TEMPERATURE (DEG F) : 56.0 

SEARCH FOR FREQUENCY PEAKS WITHIN RANGE OF 50.0 TO 250.0 HERTZ 
MICROPHONE COORDINATES : 


MIC 

# 

1 

X 

-4.0 

Y 

0.0 

MIC 

# 

2 

X 

4.0 

Y 

0.0 

MIC 

# 

3 

X 

0.0 

Y 

-6.9 

MIC 

# 

4 

X 

-2.0 

Y 

-3.4 

MIC 

# 

5 

X 

0.0 

Y 

0.0 

MIC 

# 

6 

X 

2.0 

Y 

-3.4 

MIC 

# 

7 

X 

0.0 

Y 

-3.4 

MIC 

# 

8 

X 

-1.0 

Y 

-1.7 

MIC 

# 

9 

X 

1.0 

Y 

-1.7 


***** TIFT DATA PROCESSING ***** 

SERIAL NUMBER : 15 

NUMBER OF CHANNELS : 10 

HEADER : A/C TRACKING/LANGLEY TEST NO. AC RUN NO. 5 


CHANNEL 

NAME 

UNITS 

1 

TIME 

SECONDS 

2 

MIC1 

DYNES/CM2 

3 

MIC2 

DYNES/CM2 

4 

MIC3 

DYNES/CM2 

5 

MIC4 

DYNES/CM2 

6 

MIC5 

DYNES/CM2 

7 

MIC6 

DYNES/CM2 

8 

MIC7 

DYNES/CM2 

9 

MIC8 

DYNES/CM2 

10 

MIC9 

DYNES/CM2 


START TIME = 20999.6812 


END TIME = 21000.1807 


Figure 3a: Output File BEARING 
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***** PEAK SEARCH RESULTS ***** 


PEAK # FREQUENCY (HZ) SOUND PRESSURE LEVEL (DB) 


140.0 

56.5 

116.0 

56.1 

93.0 

54.1 

163.0 

51.9 

209.0 

51.3 

186.0 

50.7 

233.0 

49.7 

70.0 

47.5 


***** SOURCE LOCATION RESULTS ***** 


FREQUENCY (HZ) # PEAKS 


140.0 2 

116.0 3 

93.0 1 

163.0 1 

209.0 3 

186.0 2 

233.0 3 

70.0 1 


AZIMUTH (DEG) 

POWER 

121.0 

56.4 

300.0 

52.0 

122.0 

56.0 

341.0 

48.9 

244.0 

48.9 

121.0 

54.1 

120.0 

51.9 

121.0 

51.2 

347.0 

47.5 

254.0 

47.3 

118.0 

50.9 

251.0 

45.7 

120.0 

49.5 

254.0 

46.6 

348.0 

46.6 

119.0 

47.5 


Figure 3b: Output File BEARING (cont.) 
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CONCLUDING REMARKS 


This paper has presented the mathematical derivation of a 
new acoustic source bearing estimation method, the ASAT 
algorithm, and described the acoustic array processing computer 
program, ASBE, which was developed to use it. By i ncorporati ng 
the high frequency resolution of the Fast Fourier Transform (FFT) 
and the bearing estimation capability of Capon's Maximum 
Likelihood Method (MLM), ASAT will estimate the bearing of an 
acoustic source while minimizing the effects of foldover aliasing 
and spatial aliasing. Unlike the MLM, ASAT does not require the 
inverse of the sensor array's cross power spectral density 
matrix. This feature allows the analysis of data which may yield 
a singular matrix. 

Using the ASBE program, a benchmark test case consisting of 
60 runs was analyzed and the results show a high accuracy in 
locating the sound source. Polar bearing plots from the first 
run show how sidelobes became more apparent when a subset of the 
sensor data was used in the analysis. Since the highest 
frequency to be analyzed is usually not known, this sidelobe 
effect is an acceptable alternative to deleting those frequencies 
which create spatial aliasing. 
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